function  [Ns2D,Ns3D,gradNs3D] = basisNj(deg)


switch (deg)
    case(1)
 
        baryc = ...
            [1 1 0 0 0;...
             2 0 1 0 0;...
             3 0 0 1 0;...
             4 0 0 0 1];
        Ns3D{1} = @(l1,l2,l3,l4) l1;
        Ns3D{2} = @(l1,l2,l3,l4) l2;
        Ns3D{3} = @(l1,l2,l3,l4) l3;
        Ns3D{4} =  @(l1,l2,l3,l4) l4;

        Ns2D{1} = @(l1,l2,l3) l1;
        Ns2D{2} = @(l1,l2,l3) l2;
        Ns2D{3} = @(l1,l2,l3) l3;
        
        
        gradNs3D = {...
            @(l1,l2,l3,l4)-1+0.*l1,...
            @(l1,l2,l3,l4)-1+0.*l1,...
            @(l1,l2,l3,l4)-1+0.*l1;...
            @(l1,l2,l3,l4)1+0.*l1,...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)0+0.*l1;...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)1+0.*l1,...
            @(l1,l2,l3,l4)0+0.*l1;...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)1+0.*l1;...
        }';

    case(2)
 
        baryc = ...
            [1        0         0        0;...
             0        1         0        0;...
             0        0         1        0;...
             0        0         0        1;...
           1/2      1/2         0        0;...
             0      1/2       1/2        0;...
           1/2        0       1/2        0;...
           1/2        0         0      1/2;...
             0        0       1/2      1/2;...
             0      1/2         0      1/2];

        Ns3D{1}  = @(l1,l2,l3,l4) 2*l1.*(l1 - 1/2);
        Ns3D{2}  = @(l1,l2,l3,l4) 2*l2.*(l2 - 1/2);
        Ns3D{3}  = @(l1,l2,l3,l4) 2*l3.*(l3 - 1/2);
        Ns3D{4}  = @(l1,l2,l3,l4) 2*l4.*(l4 - 1/2);
        Ns3D{5}  = @(l1,l2,l3,l4) 4*l1.*l2;
        Ns3D{6}  = @(l1,l2,l3,l4) 4*l2.*l3;
        Ns3D{7}  = @(l1,l2,l3,l4) 4*l1.*l3;
        Ns3D{8}  = @(l1,l2,l3,l4) 4*l1.*l4;
        Ns3D{9}  = @(l1,l2,l3,l4) 4*l3.*l4;
        Ns3D{10} = @(l1,l2,l3,l4) 4*l2.*l4;


        Ns2D{1} = @(l1,l2,l3) 2*l1.*(l1 - 1/2);
        Ns2D{2} = @(l1,l2,l3) 2*l2.*(l2 - 1/2);
        Ns2D{3} = @(l1,l2,l3) 2*l3.*(l3 - 1/2);
        Ns2D{4} = @(l1,l2,l3) 4*l1.*l2;
        Ns2D{5} = @(l1,l2,l3) 4*l2.*l3;
        Ns2D{6} = @(l1,l2,l3) 4*l1.*l3;
        
        gradNs3D = {...
            @(l1,l2,l3,l4)1 - 4.*l1,...
            @(l1,l2,l3,l4)1 - 4.*l1,...
            @(l1,l2,l3,l4)1 - 4.*l1;...
            @(l1,l2,l3,l4)4.*l2 - 1,...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)0+0.*l1;...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)4.*l3 - 1,...
            @(l1,l2,l3,l4)0+0.*l1;...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)4.*l4 - 1;...
            @(l1,l2,l3,l4)4.*l1 - 4.*l2,...
            @(l1,l2,l3,l4)-4.*l2,...
            @(l1,l2,l3,l4)-4.*l2;...
            @(l1,l2,l3,l4)4.*l3,...
            @(l1,l2,l3,l4)4.*l2,...
            @(l1,l2,l3,l4)0+0.*l1;...
            @(l1,l2,l3,l4)-4.*l3,...
            @(l1,l2,l3,l4)4.*l1 - 4.*l3,...
            @(l1,l2,l3,l4)-4.*l3;...
            @(l1,l2,l3,l4)-4.*l4,...
            @(l1,l2,l3,l4)-4.*l4,...
            @(l1,l2,l3,l4)4.*l1 - 4.*l4;...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)4.*l4,...
            @(l1,l2,l3,l4)4.*l3;...
            @(l1,l2,l3,l4)4.*l4,...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)4.*l2;...
        }';

    case(3) 
        % P3

        baryc = ...
            [1   1        0         0        0;...
            2   0        1         0        0;...
            3   0        0         1        0;...
            4   0        0         0        1;...
            %  Edges
            %  y = 0, z = 0;
            5   2/3      1/3       0        0;...
            6   1/3      2/3       0        0;...
            % z = 0, x+y =1;
            7   0        2/3      1/3       0;...
            8   0        1/3      2/3       0;...
            % x = 0,z  =0;
            9   1/3      0        2/3       0;...
            10  2/3      0        1/3       0;...
            %    x = 0, y  =0;
            11  1/3      0         0        2/3;...
            12  2/3      0         0        1/3;...
            %    x = 0, y+z  =0;
            13  0        0        1/3       2/3;...
            14  0        0         2/3      1/3;...
            %    y = 0, x+z  =1;
            15  0        1/3       0        2/3;...
            16  0        2/3       0        1/3;...
            %
            % internal faces nodes
            % z=9
            17  1/3        1/3    1/3           0;...
            % y=0
            18  1/3       1/3       0          1/3;...
            % x=0
            19  1/3         0     1/3          1/3;...
            % x+y+z=1
            20  0         1/3     1/3          1/3;...
            ];
        baryc=baryc(:,2:end);
        Ns3D{01}=@(l1,l2,l3,l4) (9*l1.*(l1 - 1/3).*(l1 - 2/3))/2;
        Ns3D{02}=@(l1,l2,l3,l4) (9*l2.*(l2 - 1/3).*(l2 - 2/3))/2;
        Ns3D{03}=@(l1,l2,l3,l4) (9*l3.*(l3 - 1/3).*(l3 - 2/3))/2;
        Ns3D{04}=@(l1,l2,l3,l4) (9*l4.*(l4 - 1/3).*(l4 - 2/3))/2;
        Ns3D{05}=@(l1,l2,l3,l4)       (27*l1.*l2.*(l1 - 1/3))/2;
        Ns3D{06}=@(l1,l2,l3,l4)       (27*l1.*l2.*(l2 - 1/3))/2;
        Ns3D{07}=@(l1,l2,l3,l4)       (27*l2.*l3.*(l2 - 1/3))/2;
        Ns3D{08}=@(l1,l2,l3,l4)       (27*l2.*l3.*(l3 - 1/3))/2;
        Ns3D{09}=@(l1,l2,l3,l4)       (27*l1.*l3.*(l3 - 1/3))/2;
        Ns3D{10}=@(l1,l2,l3,l4)       (27*l1.*l3.*(l1 - 1/3))/2;
        Ns3D{11}=@(l1,l2,l3,l4)       (27*l1.*l4.*(l4 - 1/3))/2;
        Ns3D{12}=@(l1,l2,l3,l4)       (27*l1.*l4.*(l1 - 1/3))/2;
        Ns3D{13}=@(l1,l2,l3,l4)       (27*l3.*l4.*(l4 - 1/3))/2;
        Ns3D{14}=@(l1,l2,l3,l4)       (27*l3.*l4.*(l3 - 1/3))/2;
        Ns3D{15}=@(l1,l2,l3,l4)       (27*l2.*l4.*(l4 - 1/3))/2;
        Ns3D{16}=@(l1,l2,l3,l4)       (27*l2.*l4.*(l2 - 1/3))/2;
        Ns3D{17}=@(l1,l2,l3,l4)                   27*l1.*l2.*l3;
        Ns3D{18}=@(l1,l2,l3,l4)                   27*l1.*l2.*l4;
        Ns3D{19}=@(l1,l2,l3,l4)                   27*l1.*l3.*l4;
        Ns3D{20}=@(l1,l2,l3,l4)                   27*l2.*l3.*l4;



        Ns2D{01}=@(l1,l2,l3) (9*l1.*(l1 - 1/3).*(l1 - 2/3))/2;
        Ns2D{02}=@(l1,l2,l3) (9*l2.*(l2 - 1/3).*(l2 - 2/3))/2;
        Ns2D{03}=@(l1,l2,l3) (9*l3.*(l3 - 1/3).*(l3 - 2/3))/2;
        Ns2D{04}=@(l1,l2,l3)       (27*l1.*l2.*(l1 - 1/3))/2;
        Ns2D{05}=@(l1,l2,l3)       (27*l1.*l2.*(l2 - 1/3))/2;
        Ns2D{06}=@(l1,l2,l3)       (27*l2.*l3.*(l2 - 1/3))/2;
        Ns2D{07}=@(l1,l2,l3)       (27*l2.*l3.*(l3 - 1/3))/2;
        Ns2D{08}=@(l1,l2,l3)       (27*l1.*l3.*(l3 - 1/3))/2;
        Ns2D{09}=@(l1,l2,l3)       (27*l1.*l3.*(l1 - 1/3))/2;
        Ns2D{10}=@(l1,l2,l3)                   27*l1.*l2.*l3;
        
        
        
        gradNs3D = {...
            @(l1,l2,l3,l4)9.*l1 - (27.*l1.^2)/2 - 1,...
            @(l1,l2,l3,l4)9.*l1 - (27.*l1.^2)/2 - 1,...
            @(l1,l2,l3,l4)9.*l1 - (27.*l1.^2)/2 - 1;...
            @(l1,l2,l3,l4)(27.*l2.^2)/2 - 9.*l2 + 1,...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)0+0.*l1;...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)(27.*l3.^2)/2 - 9.*l3 + 1,...
            @(l1,l2,l3,l4)0+0.*l1;...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)(27.*l4.^2)/2 - 9.*l4 + 1;...
            @(l1,l2,l3,l4)(9.*l2)/2 - (9.*l1)/2 - 27.*l1.*l2 + (27.*l1.^2)/2,...
            @(l1,l2,l3,l4)-l2.*(27.*l1 - 9/2),...
            @(l1,l2,l3,l4)-l2.*(27.*l1 - 9/2);...
            @(l1,l2,l3,l4)(9.*l2)/2 - (9.*l1)/2 + 27.*l1.*l2 - (27.*l2.^2)/2,...
            @(l1,l2,l3,l4)-(9.*l2.*(3.*l2 - 1))/2,...
            @(l1,l2,l3,l4)-(9.*l2.*(3.*l2 - 1))/2;...
            @(l1,l2,l3,l4)l3.*(27.*l2 - 9/2),...
            @(l1,l2,l3,l4)(9.*l2.*(3.*l2 - 1))/2,...
            @(l1,l2,l3,l4)0+0.*l1;...
            @(l1,l2,l3,l4)(9.*l3.*(3.*l3 - 1))/2,...
            @(l1,l2,l3,l4)(9.*l2.*(6.*l3 - 1))/2,...
            @(l1,l2,l3,l4)0+0.*l1;...
            @(l1,l2,l3,l4)-(9.*l3.*(3.*l3 - 1))/2,...
            @(l1,l2,l3,l4)(9.*l3)/2 - (9.*l1)/2 + 27.*l1.*l3 - (27.*l3.^2)/2,...
            @(l1,l2,l3,l4)-(9.*l3.*(3.*l3 - 1))/2;...
            @(l1,l2,l3,l4)-l3.*(27.*l1 - 9/2),...
            @(l1,l2,l3,l4)(9.*l3)/2 - (9.*l1)/2 - 27.*l1.*l3 + (27.*l1.^2)/2,...
            @(l1,l2,l3,l4)-l3.*(27.*l1 - 9/2);...
            @(l1,l2,l3,l4)-(9.*l4.*(3.*l4 - 1))/2,...
            @(l1,l2,l3,l4)-(9.*l4.*(3.*l4 - 1))/2,...
            @(l1,l2,l3,l4)(9.*l4)/2 - (9.*l1)/2 + 27.*l1.*l4 - (27.*l4.^2)/2;...
            @(l1,l2,l3,l4)-l4.*(27.*l1 - 9/2),...
            @(l1,l2,l3,l4)-l4.*(27.*l1 - 9/2),...
            @(l1,l2,l3,l4)(9.*l4)/2 - (9.*l1)/2 - 27.*l1.*l4 + (27.*l1.^2)/2;...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)(9.*l4.*(3.*l4 - 1))/2,...
            @(l1,l2,l3,l4)(9.*l3.*(6.*l4 - 1))/2;...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)l4.*(27.*l3 - 9/2),...
            @(l1,l2,l3,l4)(9.*l3.*(3.*l3 - 1))/2;...
            @(l1,l2,l3,l4)(9.*l4.*(3.*l4 - 1))/2,...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)(9.*l2.*(6.*l4 - 1))/2;...
            @(l1,l2,l3,l4)l4.*(27.*l2 - 9/2),...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)(9.*l2.*(3.*l2 - 1))/2;...
            @(l1,l2,l3,l4)l3.*(27.*l1 - 27.*l2),...
            @(l1,l2,l3,l4)27.*l2.*(l1 - l3),...
            @(l1,l2,l3,l4)-27.*l2.*l3;...
            @(l1,l2,l3,l4)l4.*(27.*l1 - 27.*l2),...
            @(l1,l2,l3,l4)-27.*l2.*l4,...
            @(l1,l2,l3,l4)27.*l2.*(l1 - l4);...
            @(l1,l2,l3,l4)-27.*l3.*l4,...
            @(l1,l2,l3,l4)l4.*(27.*l1 - 27.*l3),...
            @(l1,l2,l3,l4)27.*l3.*(l1 - l4);...
            @(l1,l2,l3,l4)27.*l3.*l4,...
            @(l1,l2,l3,l4)27.*l2.*l4,...
            @(l1,l2,l3,l4)27.*l2.*l3;...
        }';

      
    case(4)
        % P4 

        baryc= ...
            [1   1        0         0        0;...
             2   0        1         0        0;...
             3   0        0         1        0;...
             4   0        0         0        1;...
            %  edges
            %  y = z = 0
             5   3/4      1/4       0        0;...
             6   2/4      2/4       0        0;...
             7   1/4      3/4       0        0;...
            %    z = 0, x+y=1
             8   0        3/4      1/4       0;...
             9   0        2/4      2/4       0;...
             10  0        1/4      3/4       0;...
            %        x = z =0
             11   1/4       0      3/4       0;...
             12   2/4       0      2/4       0;...
             13   3/4       0      1/4       0;...
            %        x = 0, y  =0;
            14    1/4       0        0     3/4;...
            15    2/4       0        0     2/4;...
            16    3/4       0        0     1/4;...
            %        x= 0, y+z=0
            17  0        0       1/4        3/4;...
            18  0        0       2/4        2/4;...
            19  0        0       3/4        1/4;...
            %         y=0, x+z=1
            20  0       1/4        0      3/4;...
            21  0       2/4        0       2/4;...
            22  0       3/4        0       1 /4;...
            % internal faces z=0
            23  2/4      1/4      1/4        0;...
            24  1/4      1/4      2/4        0;
            25  1/4      2/4      1/4        0; ...
            % internal faces y =0
            26  2/4      1/4       0     1/4     ;...
            27  1/4      2/4       0     1/4     ;...
            28  1/4      1/4       0     2/4     ;...
            % internal  faces x= 0
            29           2/4 0    1/4       1/4       ;...
            30           1/4 0     1/4      2/4    ;...
            31           1/4 0     2/4      1/4    ;...
            % internal faces x+y+z=1
            32  0      1/4        1/4      2/4    ;...
            33  0      2/4        1/4      1/4    ;...
            34  0      1/4        2/4      1/4    ;...
            % internal  nodes
            35  1/4      1/4     1/4       1/4;...
            ];

        baryc = baryc(:,2:end);
        Ns3D{01}= @ (l1,l2,l3,l4) (32*l1.*(l1 - 1/2).*(l1 - 1/4).*(l1 - 3/4))/3;
        Ns3D{02}= @ (l1,l2,l3,l4) (32*l2.*(l2 - 1/2).*(l2 - 1/4).*(l2 - 3/4))/3;
        Ns3D{03}= @ (l1,l2,l3,l4) (32*l3.*(l3 - 1/2).*(l3 - 1/4).*(l3 - 3/4))/3;
        Ns3D{04}= @ (l1,l2,l3,l4) (32*l4.*(l4 - 1/2).*(l4 - 1/4).*(l4 - 3/4))/3;
        Ns3D{05}= @ (l1,l2,l3,l4)        (128*l1.*l2.*(l1 - 1/2).*(l1 - 1/4))/3;
        Ns3D{06}= @ (l1,l2,l3,l4)             64*l1.*l2.*(l1 - 1/4).*(l2 - 1/4);
        Ns3D{07}= @ (l1,l2,l3,l4)        (128*l1.*l2.*(l2 - 1/2).*(l2 - 1/4))/3;
        Ns3D{08}= @ (l1,l2,l3,l4)        (128*l2.*l3.*(l2 - 1/2).*(l2 - 1/4))/3;
        Ns3D{09}= @ (l1,l2,l3,l4)             64*l2.*l3.*(l2 - 1/4).*(l3 - 1/4);
        Ns3D{10}= @ (l1,l2,l3,l4)        (128*l2.*l3.*(l3 - 1/2).*(l3 - 1/4))/3;
        Ns3D{11}= @ (l1,l2,l3,l4)       (128*l1.*l3.*(l3 - 1/2).*(l3 - 1/4))/3 ;
        Ns3D{12}= @ (l1,l2,l3,l4)            64*l1.*l3.*(l1 - 1/4).*(l3 - 1/4) ;
        Ns3D{13}= @ (l1,l2,l3,l4)       (128*l1.*l3.*(l1 - 1/2).*(l1 - 1/4))/3 ;
        Ns3D{14}= @ (l1,l2,l3,l4)       (128*l1.*l4.*(l4 - 1/2).*(l4 - 1/4))/3 ;
        Ns3D{15}= @ (l1,l2,l3,l4)            64*l1.*l4.*(l1 - 1/4).*(l4 - 1/4) ;
        Ns3D{16}= @ (l1,l2,l3,l4)       (128*l1.*l4.*(l1 - 1/2).*(l1 - 1/4))/3 ;
        Ns3D{17}= @ (l1,l2,l3,l4)       (128*l3.*l4.*(l4 - 1/2).*(l4 - 1/4))/3 ;
        Ns3D{18}= @ (l1,l2,l3,l4)            64*l3.*l4.*(l3 - 1/4).*(l4 - 1/4) ;
        Ns3D{19}= @ (l1,l2,l3,l4)       (128*l3.*l4.*(l3 - 1/2).*(l3 - 1/4))/3 ;
        Ns3D{20}= @ (l1,l2,l3,l4)       (128*l2.*l4.*(l4 - 1/2).*(l4 - 1/4))/3 ;
        Ns3D{21}= @ (l1,l2,l3,l4)            64*l2.*l4.*(l2 - 1/4).*(l4 - 1/4) ;
        Ns3D{22}= @ (l1,l2,l3,l4)       (128*l2.*l4.*(l2 - 1/2).*(l2 - 1/4))/3 ;
        Ns3D{23}= @ (l1,l2,l3,l4)                   128*l1.*l2.*l3.*(l1 - 1/4) ;
        Ns3D{24}= @ (l1,l2,l3,l4)                   128*l1.*l2.*l3.*(l3 - 1/4) ;
        Ns3D{25}= @ (l1,l2,l3,l4)                   128*l1.*l2.*l3.*(l2 - 1/4) ;
        Ns3D{26}= @ (l1,l2,l3,l4)                   128*l1.*l2.*l4.*(l1 - 1/4) ;
        Ns3D{27}= @ (l1,l2,l3,l4)                   128*l1.*l2.*l4.*(l2 - 1/4) ;
        Ns3D{28}= @ (l1,l2,l3,l4)                   128*l1.*l2.*l4.*(l4 - 1/4) ;
        Ns3D{29}= @ (l1,l2,l3,l4)                   128*l1.*l3.*l4.*(l1 - 1/4) ;
        Ns3D{30}= @ (l1,l2,l3,l4)                   128*l1.*l3.*l4.*(l4 - 1/4) ;
        Ns3D{31}= @ (l1,l2,l3,l4)                   128*l1.*l3.*l4.*(l3 - 1/4) ;
        Ns3D{32}= @ (l1,l2,l3,l4)                   128*l2.*l3.*l4.*(l4 - 1/4) ;
        Ns3D{33}= @ (l1,l2,l3,l4)                   128*l2.*l3.*l4.*(l2 - 1/4) ;
        Ns3D{34}= @ (l1,l2,l3,l4)                   128*l2.*l3.*l4.*(l3 - 1/4) ;
        Ns3D{35}= @ (l1,l2,l3,l4)                           256*l1.*l2.*l3.*l4 ;



        Ns2D{01}= @ (l1,l2,l3) (32*l1.*(l1 - 1/2).*(l1 - 1/4).*(l1 - 3/4))/3;
        Ns2D{02}= @ (l1,l2,l3) (32*l2.*(l2 - 1/2).*(l2 - 1/4).*(l2 - 3/4))/3;
        Ns2D{03}= @ (l1,l2,l3) (32*l3.*(l3 - 1/2).*(l3 - 1/4).*(l3 - 3/4))/3;
        Ns2D{04}= @ (l1,l2,l3)        (128*l1.*l2.*(l1 - 1/2).*(l1 - 1/4))/3;
        Ns2D{05}= @ (l1,l2,l3)             64*l1.*l2.*(l1 - 1/4).*(l2 - 1/4);
        Ns2D{06}= @ (l1,l2,l3)        (128*l1.*l2.*(l2 - 1/2).*(l2 - 1/4))/3;
        Ns2D{07}= @ (l1,l2,l3)        (128*l2.*l3.*(l2 - 1/2).*(l2 - 1/4))/3;
        Ns2D{08}= @ (l1,l2,l3)             64*l2.*l3.*(l2 - 1/4).*(l3 - 1/4);
        Ns2D{09}= @ (l1,l2,l3)        (128*l2.*l3.*(l3 - 1/2).*(l3 - 1/4))/3;
        Ns2D{10}= @ (l1,l2,l3)       (128*l1.*l3.*(l3 - 1/2).*(l3 - 1/4))/3 ;
        Ns2D{11}= @ (l1,l2,l3)            64*l1.*l3.*(l1 - 1/4).*(l3 - 1/4) ;
        Ns2D{12}= @ (l1,l2,l3)       (128*l1.*l3.*(l1 - 1/2).*(l1 - 1/4))/3 ;
        Ns2D{13}= @ (l1,l2,l3)                   128*l1.*l2.*l3.*(l1 - 1/4) ;
        Ns2D{14}= @ (l1,l2,l3)                   128*l1.*l2.*l3.*(l2 - 1/4) ;
        Ns2D{15}= @ (l1,l2,l3)                   128*l1.*l2.*l3.*(l3 - 1/4) ;
        
        
        
        gradNs3D = {...
            @(l1,l2,l3,l4)-((8.*l1)/3 - 1).*(16.*l1.^2 - 12.*l1 + 1),...
            @(l1,l2,l3,l4)-((8.*l1)/3 - 1).*(16.*l1.^2 - 12.*l1 + 1),...
            @(l1,l2,l3,l4)-((8.*l1)/3 - 1).*(16.*l1.^2 - 12.*l1 + 1);...
            @(l1,l2,l3,l4)((8.*l2)/3 - 1).*(16.*l2.^2 - 12.*l2 + 1),...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)0+0.*l1;...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)((8.*l3)/3 - 1).*(16.*l3.^2 - 12.*l3 + 1),...
            @(l1,l2,l3,l4)0+0.*l1;...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)((8.*l4)/3 - 1).*(16.*l4.^2 - 12.*l4 + 1);...
            @(l1,l2,l3,l4)(16.*l1)/3 - (16.*l2)/3 + 64.*l1.*l2 - 128.*l1.^2.*l2 - 32.*l1.^2 + (128.*l1.^3)/3,...
            @(l1,l2,l3,l4)-l2.*(128.*l1.^2 - 64.*l1 + 16/3),...
            @(l1,l2,l3,l4)-l2.*(128.*l1.^2 - 64.*l1 + 16/3);...
            @(l1,l2,l3,l4)-(4.*l1 - 4.*l2).*(4.*l1 + 4.*l2 - 32.*l1.*l2 - 1),...
            @(l1,l2,l3,l4)-l2.*(4.*l2 - 1).*(32.*l1 - 4),...
            @(l1,l2,l3,l4)-l2.*(4.*l2 - 1).*(32.*l1 - 4);...
            @(l1,l2,l3,l4)(16.*l1)/3 - (16.*l2)/3 - 64.*l1.*l2 + 128.*l1.*l2.^2 + 32.*l2.^2 - (128.*l2.^3)/3,...
            @(l1,l2,l3,l4)-(16.*l2.*(2.*l2 - 1).*(4.*l2 - 1))/3,...
            @(l1,l2,l3,l4)-(16.*l2.*(2.*l2 - 1).*(4.*l2 - 1))/3;...
            @(l1,l2,l3,l4)l3.*(128.*l2.^2 - 64.*l2 + 16/3),...
            @(l1,l2,l3,l4)(16.*l2.*(2.*l2 - 1).*(4.*l2 - 1))/3,...
            @(l1,l2,l3,l4)0+0.*l1;...
            @(l1,l2,l3,l4)l3.*(4.*l3 - 1).*(32.*l2 - 4),...
            @(l1,l2,l3,l4)4.*l2.*(4.*l2 - 1).*(8.*l3 - 1),...
            @(l1,l2,l3,l4)0+0.*l1;...
            @(l1,l2,l3,l4)(16.*l3.*(2.*l3 - 1).*(4.*l3 - 1))/3,...
            @(l1,l2,l3,l4)(16.*l2.*(24.*l3.^2 - 12.*l3 + 1))/3,...
            @(l1,l2,l3,l4)0+0.*l1;...
            @(l1,l2,l3,l4)-(16.*l3.*(2.*l3 - 1).*(4.*l3 - 1))/3,...
            @(l1,l2,l3,l4)(16.*l1)/3 - (16.*l3)/3 - 64.*l1.*l3 + 128.*l1.*l3.^2 + 32.*l3.^2 - (128.*l3.^3)/3,...
            @(l1,l2,l3,l4)-(16.*l3.*(2.*l3 - 1).*(4.*l3 - 1))/3;...
            @(l1,l2,l3,l4)-l3.*(4.*l3 - 1).*(32.*l1 - 4),...
            @(l1,l2,l3,l4)-(4.*l1 - 4.*l3).*(4.*l1 + 4.*l3 - 32.*l1.*l3 - 1),...
            @(l1,l2,l3,l4)-l3.*(4.*l3 - 1).*(32.*l1 - 4);...
            @(l1,l2,l3,l4)-l3.*(128.*l1.^2 - 64.*l1 + 16/3),...
            @(l1,l2,l3,l4)(16.*l1)/3 - (16.*l3)/3 + 64.*l1.*l3 - 128.*l1.^2.*l3 - 32.*l1.^2 + (128.*l1.^3)/3,...
            @(l1,l2,l3,l4)-l3.*(128.*l1.^2 - 64.*l1 + 16/3);...
            @(l1,l2,l3,l4)-(16.*l4.*(2.*l4 - 1).*(4.*l4 - 1))/3,...
            @(l1,l2,l3,l4)-(16.*l4.*(2.*l4 - 1).*(4.*l4 - 1))/3,...
            @(l1,l2,l3,l4)(16.*l1)/3 - (16.*l4)/3 - 64.*l1.*l4 + 128.*l1.*l4.^2 + 32.*l4.^2 - (128.*l4.^3)/3;...
            @(l1,l2,l3,l4)-l4.*(4.*l4 - 1).*(32.*l1 - 4),...
            @(l1,l2,l3,l4)-l4.*(4.*l4 - 1).*(32.*l1 - 4),...
            @(l1,l2,l3,l4)-(4.*l1 - 4.*l4).*(4.*l1 + 4.*l4 - 32.*l1.*l4 - 1);...
            @(l1,l2,l3,l4)-l4.*(128.*l1.^2 - 64.*l1 + 16/3),...
            @(l1,l2,l3,l4)-l4.*(128.*l1.^2 - 64.*l1 + 16/3),...
            @(l1,l2,l3,l4)(16.*l1)/3 - (16.*l4)/3 + 64.*l1.*l4 - 128.*l1.^2.*l4 - 32.*l1.^2 + (128.*l1.^3)/3;...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)(16.*l4.*(2.*l4 - 1).*(4.*l4 - 1))/3,...
            @(l1,l2,l3,l4)(16.*l3.*(24.*l4.^2 - 12.*l4 + 1))/3;...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)l4.*(4.*l4 - 1).*(32.*l3 - 4),...
            @(l1,l2,l3,l4)4.*l3.*(4.*l3 - 1).*(8.*l4 - 1);...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)l4.*(128.*l3.^2 - 64.*l3 + 16/3),...
            @(l1,l2,l3,l4)(16.*l3.*(2.*l3 - 1).*(4.*l3 - 1))/3;...
            @(l1,l2,l3,l4)(16.*l4.*(2.*l4 - 1).*(4.*l4 - 1))/3,...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)(16.*l2.*(24.*l4.^2 - 12.*l4 + 1))/3;...
            @(l1,l2,l3,l4)l4.*(4.*l4 - 1).*(32.*l2 - 4),...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)4.*l2.*(4.*l2 - 1).*(8.*l4 - 1);...
            @(l1,l2,l3,l4)l4.*(128.*l2.^2 - 64.*l2 + 16/3),...
            @(l1,l2,l3,l4)0+0.*l1,...
            @(l1,l2,l3,l4)(16.*l2.*(2.*l2 - 1).*(4.*l2 - 1))/3;...
            @(l1,l2,l3,l4)-l3.*(32.*l1 - 32.*l2 + 256.*l1.*l2 - 128.*l1.^2),...
            @(l1,l2,l3,l4)-32.*l2.*(l1 - l3 + 8.*l1.*l3 - 4.*l1.^2),...
            @(l1,l2,l3,l4)-l2.*l3.*(256.*l1 - 32);...
            @(l1,l2,l3,l4)l3.*(4.*l3 - 1).*(32.*l1 - 32.*l2),...
            @(l1,l2,l3,l4)-32.*l2.*(l1 - l3 - 8.*l1.*l3 + 4.*l3.^2),...
            @(l1,l2,l3,l4)-32.*l2.*l3.*(4.*l3 - 1);...
            @(l1,l2,l3,l4)-l3.*(32.*l1 - 32.*l2 - 256.*l1.*l2 + 128.*l2.^2),...
            @(l1,l2,l3,l4)32.*l2.*(4.*l2 - 1).*(l1 - l3),...
            @(l1,l2,l3,l4)-32.*l2.*l3.*(4.*l2 - 1);...
            @(l1,l2,l3,l4)-l4.*(32.*l1 - 32.*l2 + 256.*l1.*l2 - 128.*l1.^2),...
            @(l1,l2,l3,l4)-l2.*l4.*(256.*l1 - 32),...
            @(l1,l2,l3,l4)-32.*l2.*(l1 - l4 + 8.*l1.*l4 - 4.*l1.^2);...
            @(l1,l2,l3,l4)-l4.*(32.*l1 - 32.*l2 - 256.*l1.*l2 + 128.*l2.^2),...
            @(l1,l2,l3,l4)-32.*l2.*l4.*(4.*l2 - 1),...
            @(l1,l2,l3,l4)32.*l2.*(4.*l2 - 1).*(l1 - l4);...
            @(l1,l2,l3,l4)l4.*(4.*l4 - 1).*(32.*l1 - 32.*l2),...
            @(l1,l2,l3,l4)-32.*l2.*l4.*(4.*l4 - 1),...
            @(l1,l2,l3,l4)-32.*l2.*(l1 - l4 - 8.*l1.*l4 + 4.*l4.^2);...
            @(l1,l2,l3,l4)-l3.*l4.*(256.*l1 - 32),...
            @(l1,l2,l3,l4)-l4.*(32.*l1 - 32.*l3 + 256.*l1.*l3 - 128.*l1.^2),...
            @(l1,l2,l3,l4)-32.*l3.*(l1 - l4 + 8.*l1.*l4 - 4.*l1.^2);...
            @(l1,l2,l3,l4)-32.*l3.*l4.*(4.*l4 - 1),...
            @(l1,l2,l3,l4)l4.*(4.*l4 - 1).*(32.*l1 - 32.*l3),...
            @(l1,l2,l3,l4)-32.*l3.*(l1 - l4 - 8.*l1.*l4 + 4.*l4.^2);...
            @(l1,l2,l3,l4)-32.*l3.*l4.*(4.*l3 - 1),...
            @(l1,l2,l3,l4)-l4.*(32.*l1 - 32.*l3 - 256.*l1.*l3 + 128.*l3.^2),...
            @(l1,l2,l3,l4)32.*l3.*(4.*l3 - 1).*(l1 - l4);...
            @(l1,l2,l3,l4)32.*l3.*l4.*(4.*l4 - 1),...
            @(l1,l2,l3,l4)32.*l2.*l4.*(4.*l4 - 1),...
            @(l1,l2,l3,l4)32.*l2.*l3.*(8.*l4 - 1);...
            @(l1,l2,l3,l4)l3.*l4.*(256.*l2 - 32),...
            @(l1,l2,l3,l4)32.*l2.*l4.*(4.*l2 - 1),...
            @(l1,l2,l3,l4)32.*l2.*l3.*(4.*l2 - 1);...
            @(l1,l2,l3,l4)32.*l3.*l4.*(4.*l3 - 1),...
            @(l1,l2,l3,l4)32.*l2.*l4.*(8.*l3 - 1),...
            @(l1,l2,l3,l4)32.*l2.*l3.*(4.*l3 - 1);...
            @(l1,l2,l3,l4)l3.*l4.*(256.*l1 - 256.*l2),...
            @(l1,l2,l3,l4)256.*l2.*l4.*(l1 - l3),...
            @(l1,l2,l3,l4)256.*l2.*l3.*(l1 - l4);...
        }';





end

end
